/*
 * Copyright 2010-2012 Susanta Tewari. <freecode4susant@users.sourceforge.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package provider.impl.genomemap.lc;

import commons.core.LikelihoodComputer;
import commons.util.ArrayUtil;
import genomemap.data.KSA00Data;
import org.apache.commons.lang.ArrayUtils;

import java.math.BigDecimal;

/**
 * @author Susanta Tewari
 * @since Aug 21, 2011
 */
final class KSA00LC implements LikelihoodComputer<KSA00Data> {

    /**
     * Method description
     *
     *
     * @param s
     *
     * @return
     */
    @Override
    public BigDecimal computeLogLHood(final KSA00Data s) {

        final int probe_count            = s.getModel().getProbes().size();
        final double[] probeSpacings =
            ArrayUtils.toPrimitive(s.getModel().getProbeSpacing().toArray(new Double[0]));
        final int ch_length              = s.getModel().getChromosomeLength();
        final int clone_length           = s.getModel().getCloneLength();
        final double false_positive_prob = s.getModel().getFalsePositiveProbability();
        final double false_negative_prob = s.getModel().getFalseNegativeProbability();
        final int clone_count            = s.getSampleSize();
        int counter                      = 0;
        final byte[][] data              = new byte[s.getSampleSize()][];

        for (final KSA00Data.Hybridization val : s) {
            data[counter++] = val.getScore();
        }


        // computation starts here.
        final int[] probeOrder = ArrayUtil.createSequence(s.getModel().getProbes().size(), false);

        /**
         * aa is represented by matrix a(i,j) in the paper.
         * DEPENDENT PARAMS: probe_count,clone_count,data,probe_false_pos,probe_false_neg
         */
        final double[][] aa = new double[clone_count][probe_count + 1];

        for (int cloneIndex = 0; cloneIndex < clone_count; cloneIndex++) {

            for (int probeIndex = 0; probeIndex < probe_count + 1; probeIndex++) {

                if (probeIndex == 0) {
                    aa[cloneIndex][probeIndex] = 0.0;
                } else if (data[cloneIndex][probeIndex - 1] == 0) {
                    aa[cloneIndex][probeIndex] = false_negative_prob / (1 - false_positive_prob);
                } else {
                    aa[cloneIndex][probeIndex] = (1 - false_negative_prob) / false_positive_prob;
                }
            }
        }

        /**
         * P = total # of 1`s in data: used in calculation
         * DEPENDENT PARAMS: probe_count, clone_count, data
         */
        int P = 0;

        for (int cloneIndex = 0; cloneIndex < clone_count; cloneIndex++) {

            for (int probeIndex = 0; probeIndex < probe_count; probeIndex++) {

                if (data[cloneIndex][probeIndex] == 1) {
                    P++;
                }
            }
        }

        /**
         * gap = total gap not filled by the probes; const1: used in the calculation
         * DEPENDENT PARAMS: ch_length, clone_length, probe_count, clone_count, probe_false_pos
         */
        final double const1 = clone_count * Math.log(ch_length - clone_length)
                              - P * Math.log(false_positive_prob / (1 - false_positive_prob))
                              - probe_count * clone_count * Math.log(1 - false_positive_prob);

        /**
         * R = vector of constants used in the calculation
         * DEPENDENT PARAMS: ch_length,clone_length,probe_count,clone_count,aa,probeOrder
         */
        final double[] R = PLHoodHelper.getR(ch_length, clone_length, probe_count, clone_count, aa,
                               probeOrder);
        final double val = -PLHoodHelper.Conj(probeSpacings, probe_count, clone_count, aa,
                               probeOrder, clone_length, const1, R);

        return new BigDecimal(val);
    }
}
